
setwd("INSERT WORKING DIRECTORY HERE")
library(arm)
library(simcf) #obtain from http://faculty.washington.edu/cadolph/index.php?page=60
library(tile) #obtain from http://faculty.washington.edu/cadolph/index.php?page=60


data<-read.csv("finaldataISQ.csv", header=TRUE)

#SPLIT SAMPLES

#1. CORR08 split, country fe
data_hi1<-data[data$CORR08>=median(tapply(data$CORR08, data$COUNTRY, mean), na.rm=T),]
data_lo1<-data[data$CORR08<median(tapply(data$CORR08, data$COUNTRY, mean), na.rm=T),]

formula_hi1 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi1[,25:54][,apply(data_hi1[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo1 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo1[,25:54][,apply(data_lo1[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi1<-glm(formula_hi1, data=data_hi1, family=binomial(link=logit))
model_lo1<-glm(formula_lo1, data=data_lo1, family=binomial(link=logit))

#2. CPI split, country fe
data_hi2<-data[data$CPI08>=median(tapply(data$CPI08, data$COUNTRY, mean), na.rm=T),]
data_lo2<-data[data$CPI08<median(tapply(data$CPI08, data$COUNTRY, mean), na.rm=T),]

formula_hi2 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi2[,25:54][,apply(data_hi2[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo2 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo2[,25:54][,apply(data_lo2[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi2<-glm(formula_hi2, data=data_hi2, family=binomial(link=logit))
model_lo2<-glm(formula_lo2, data=data_lo2, family=binomial(link=logit))

#3. GOVEFF split, country fe
data_hi3<-data[data$GOVEFF08>=median(tapply(data$GOVEFF08, data$COUNTRY, mean), na.rm=T),]
data_lo3<-data[data$GOVEFF08<median(tapply(data$GOVEFF08, data$COUNTRY, mean), na.rm=T),]

formula_hi3 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi3[,25:54][,apply(data_hi3[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo3 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo3[,25:54][,apply(data_lo3[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi3<-glm(formula_hi3, data=data_hi3, family=binomial(link=logit))
model_lo3<-glm(formula_lo3, data=data_lo3, family=binomial(link=logit))

#4. NIT split, country fe
data_hi4<-data[data$NITCORR>=median(tapply(data$NITCORR, data$COUNTRY, mean), na.rm=T),]
data_lo4<-data[data$NITCORR<median(tapply(data$NITCORR, data$COUNTRY, mean), na.rm=T),]

formula_hi4 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi4[,25:54][,apply(data_hi4[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo4 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo4[,25:54][,apply(data_lo4[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi4<-glm(formula_hi4, data=data_hi4, family=binomial(link=logit))
model_lo4<-glm(formula_lo4, data=data_lo4, family=binomial(link=logit))



#HIERARCHICAL MODELS

model_all0<-glmer(ISO_dunno0~ -1 + own_foreign+sales_export +logTIMEREG+
how_old+firmsize+subsidiary+stateowned+SERVICE+CORR08+
GDPCAP+DEMOC+SOVIET+EU+FDI_GDP+EXP_GDP+
(1 | COUNTRY), data=data, family=binomial(link=logit), REML=FALSE)

model_all1<-glmer(ISO_dunno0~ -1 + own_foreign+sales_export +logTIMEREG+
how_old+firmsize+subsidiary+stateowned+SERVICE+CORR08 + own_foreign*CORR08+
GDPCAP+DEMOC+SOVIET+EU+FDI_GDP+EXP_GDP+
(1 | COUNTRY), data=data, family=binomial(link=logit), REML=FALSE)

model_all2<-glmer(ISO_dunno0~ -1 + own_foreign+sales_export +logTIMEREG+
how_old+firmsize+subsidiary+stateowned+SERVICE+CORR08 + sales_export*CORR08+
GDPCAP+DEMOC+SOVIET+EU+FDI_GDP+EXP_GDP+
(1 | COUNTRY), data=data, family=binomial(link=logit), REML=FALSE)

model_all3<-glmer(ISO_dunno0~ -1 + own_foreign+sales_export +logTIMEREG+
how_old+firmsize+subsidiary+stateowned+SERVICE+CORR08 + logTIMEREG*CORR08+
GDPCAP+DEMOC+SOVIET+EU+FDI_GDP+EXP_GDP+
(1 | COUNTRY), data=data, family=binomial(link=logit), REML=FALSE)





#ROBUSTNESS CHECKS
#OMIT EU
data_noEU<-data[data$EU==0,]

#1. CORR08 split, country fe
data_noEU_hi1<-data_noEU[data_noEU$CORR08>=median(tapply(data_noEU$CORR08, data_noEU$COUNTRY, mean), na.rm=T),]
data_noEU_lo1<-data_noEU[data_noEU$CORR08<median(tapply(data_noEU$CORR08, data_noEU$COUNTRY, mean), na.rm=T),]

formula_hi1 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_noEU_hi1[,25:54][,apply(data_noEU_hi1[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo1 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_noEU_lo1[,25:54][,apply(data_noEU_lo1[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi1<-glm(formula_hi1, data=data_noEU_hi1, family=binomial(link=logit))
model_lo1<-glm(formula_lo1, data=data_noEU_lo1, family=binomial(link=logit))

#2. CPI split, country fe
data_noEU_hi2<-data_noEU[data_noEU$CPI08>=median(tapply(data_noEU$CPI08, data_noEU$COUNTRY, mean), na.rm=T),]
data_noEU_lo2<-data_noEU[data_noEU$CPI08<median(tapply(data_noEU$CPI08, data_noEU$COUNTRY, mean), na.rm=T),]

formula_hi2 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export  +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_noEU_hi2[,25:54][,apply(data_noEU_hi2[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo2 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export  +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_noEU_lo2[,25:54][,apply(data_noEU_lo2[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi2<-glm(formula_hi2, data=data_noEU_hi2, family=binomial(link=logit))
model_lo2<-glm(formula_lo2, data=data_noEU_lo2, family=binomial(link=logit))

#3. GOVEFF split, country fe
data_noEU_hi3<-data_noEU[data_noEU$GOVEFF08>=median(tapply(data_noEU$GOVEFF08, data_noEU$COUNTRY, mean), na.rm=T),]
data_noEU_lo3<-data_noEU[data_noEU$GOVEFF08<median(tapply(data_noEU$GOVEFF08, data_noEU$COUNTRY, mean), na.rm=T),]

formula_hi3 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export  +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_noEU_hi3[,25:54][,apply(data_noEU_hi3[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo3 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export  +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_noEU_lo3[,25:54][,apply(data_noEU_lo3[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi3<-glm(formula_hi3, data=data_noEU_hi3, family=binomial(link=logit))
model_lo3<-glm(formula_lo3, data=data_noEU_lo3, family=binomial(link=logit))

#4. NIT split, country fe
data_noEU_hi4<-data_noEU[data_noEU$NITCORR>=median(tapply(data_noEU$NITCORR, data_noEU$COUNTRY, mean), na.rm=T),]
data_noEU_lo4<-data_noEU[data_noEU$NITCORR<median(tapply(data_noEU$NITCORR, data_noEU$COUNTRY, mean), na.rm=T),]

formula_hi4 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export  +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_noEU_hi4[,25:54][,apply(data_noEU_hi4[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo4 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export  +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_noEU_lo4[,25:54][,apply(data_noEU_lo4[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi4<-glm(formula_hi4, data=data_noEU_hi4, family=binomial(link=logit))
model_lo4<-glm(formula_lo4, data=data_noEU_lo4, family=binomial(link=logit))



#USE ISO: PROCNO
#(counts "in process" as no cert)
#1. CORR08 split, country fe
data_hi1<-data[data$CORR08>=median(tapply(data$CORR08, data$COUNTRY, mean), na.rm=T),]
data_lo1<-data[data$CORR08<median(tapply(data$CORR08, data$COUNTRY, mean), na.rm=T),]

formula_hi1 <- paste("ISO_procNO~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi1[,25:54][,apply(data_hi1[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo1 <- paste("ISO_procNO~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo1[,25:54][,apply(data_lo1[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi1<-glm(formula_hi1, data=data_hi1, family=binomial(link=logit))
model_lo1<-glm(formula_lo1, data=data_lo1, family=binomial(link=logit))

#2. CPI split, country fe
data_hi2<-data[data$CPI08>=median(tapply(data$CPI08, data$COUNTRY, mean), na.rm=T),]
data_lo2<-data[data$CPI08<median(tapply(data$CPI08, data$COUNTRY, mean), na.rm=T),]

formula_hi2 <- paste("ISO_procNO~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi2[,25:54][,apply(data_hi2[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo2 <- paste("ISO_procNO~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo2[,25:54][,apply(data_lo2[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi2<-glm(formula_hi2, data=data_hi2, family=binomial(link=logit))
model_lo2<-glm(formula_lo2, data=data_lo2, family=binomial(link=logit))

#3. GOVEFF split, country fe
data_hi3<-data[data$GOVEFF08>=median(tapply(data$GOVEFF08, data$COUNTRY, mean), na.rm=T),]
data_lo3<-data[data$GOVEFF08<median(tapply(data$GOVEFF08, data$COUNTRY, mean), na.rm=T),]

formula_hi3 <- paste("ISO_procNO~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi3[,25:54][,apply(data_hi3[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo3 <- paste("ISO_procNO~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo3[,25:54][,apply(data_lo3[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi3<-glm(formula_hi3, data=data_hi3, family=binomial(link=logit))
model_lo3<-glm(formula_lo3, data=data_lo3, family=binomial(link=logit))

#4. NIT split, country fe
data_hi4<-data[data$NITCORR>=median(tapply(data$NITCORR, data$COUNTRY, mean), na.rm=T),]
data_lo4<-data[data$NITCORR<median(tapply(data$NITCORR, data$COUNTRY, mean), na.rm=T),]

formula_hi4 <- paste("ISO_procNO~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi4[,25:54][,apply(data_hi4[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo4 <- paste("ISO_procNO~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo4[,25:54][,apply(data_lo4[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi4<-glm(formula_hi4, data=data_hi4, family=binomial(link=logit))
model_lo4<-glm(formula_lo4, data=data_lo4, family=binomial(link=logit))



#USE ISO: ISO_dNA
#(counts "dont know" as NA instead of 0)
#1. CORR08 split, country fe
data_hi1<-data[data$CORR08>=median(tapply(data$CORR08, data$COUNTRY, mean), na.rm=T),]
data_lo1<-data[data$CORR08<median(tapply(data$CORR08, data$COUNTRY, mean), na.rm=T),]

formula_hi1 <- paste("ISO_dNA~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi1[,25:54][,apply(data_hi1[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo1 <- paste("ISO_dNA~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo1[,25:54][,apply(data_lo1[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi1<-glm(formula_hi1, data=data_hi1, family=binomial(link=logit))
model_lo1<-glm(formula_lo1, data=data_lo1, family=binomial(link=logit))

#2. CPI split, country fe
data_hi2<-data[data$CPI08>=median(tapply(data$CPI08, data$COUNTRY, mean), na.rm=T),]
data_lo2<-data[data$CPI08<median(tapply(data$CPI08, data$COUNTRY, mean), na.rm=T),]

formula_hi2 <- paste("ISO_dNA~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi2[,25:54][,apply(data_hi2[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo2 <- paste("ISO_dNA~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo2[,25:54][,apply(data_lo2[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi2<-glm(formula_hi2, data=data_hi2, family=binomial(link=logit))
model_lo2<-glm(formula_lo2, data=data_lo2, family=binomial(link=logit))

#3. GOVEFF split, country fe
data_hi3<-data[data$GOVEFF08>=median(tapply(data$GOVEFF08, data$COUNTRY, mean), na.rm=T),]
data_lo3<-data[data$GOVEFF08<median(tapply(data$GOVEFF08, data$COUNTRY, mean), na.rm=T),]

formula_hi3 <- paste("ISO_dNA~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi3[,25:54][,apply(data_hi3[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo3 <- paste("ISO_dNA~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo3[,25:54][,apply(data_lo3[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi3<-glm(formula_hi3, data=data_hi3, family=binomial(link=logit))
model_lo3<-glm(formula_lo3, data=data_lo3, family=binomial(link=logit))

#4. NIT split, country fe
data_hi4<-data[data$NITCORR>=median(tapply(data$NITCORR, data$COUNTRY, mean), na.rm=T),]
data_lo4<-data[data$NITCORR<median(tapply(data$NITCORR, data$COUNTRY, mean), na.rm=T),]

formula_hi4 <- paste("ISO_dNA~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi4[,25:54][,apply(data_hi4[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo4 <- paste("ISO_dNA~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo4[,25:54][,apply(data_lo4[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi4<-glm(formula_hi4, data=data_hi4, family=binomial(link=logit))
model_lo4<-glm(formula_lo4, data=data_lo4, family=binomial(link=logit))



#### exports/foreign>0

#1. CORR08 split, country fe
data_hi1<-data[data$CORR08>=median(tapply(data$CORR08, data$COUNTRY, mean), na.rm=T),]
data_lo1<-data[data$CORR08<median(tapply(data$CORR08, data$COUNTRY, mean), na.rm=T),]

formula_hi1 <- paste("ISO_dunno0~", "-1 + as.numeric(own_foreign>0)+as.numeric(sales_export>0) +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi1[,25:54][,apply(data_hi1[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo1 <- paste("ISO_dunno0~", "-1 + as.numeric(own_foreign>0)+as.numeric(sales_export>0) +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo1[,25:54][,apply(data_lo1[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi1<-glm(formula_hi1, data=data_hi1, family=binomial(link=logit))
model_lo1<-glm(formula_lo1, data=data_lo1, family=binomial(link=logit))

#2. CPI split, country fe
data_hi2<-data[data$CPI08>=median(tapply(data$CPI08, data$COUNTRY, mean), na.rm=T),]
data_lo2<-data[data$CPI08<median(tapply(data$CPI08, data$COUNTRY, mean), na.rm=T),]

formula_hi2 <- paste("ISO_dunno0~", "-1 + as.numeric(own_foreign>0)+as.numeric(sales_export>0) +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi2[,25:54][,apply(data_hi2[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo2 <- paste("ISO_dunno0~", "-1 + as.numeric(own_foreign>0)+as.numeric(sales_export>0) +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo2[,25:54][,apply(data_lo2[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi2<-glm(formula_hi2, data=data_hi2, family=binomial(link=logit))
model_lo2<-glm(formula_lo2, data=data_lo2, family=binomial(link=logit))

#3. GOVEFF split, country fe
data_hi3<-data[data$GOVEFF08>=median(tapply(data$GOVEFF08, data$COUNTRY, mean), na.rm=T),]
data_lo3<-data[data$GOVEFF08<median(tapply(data$GOVEFF08, data$COUNTRY, mean), na.rm=T),]

formula_hi3 <- paste("ISO_dunno0~", "-1 + as.numeric(own_foreign>0)+as.numeric(sales_export>0) +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi3[,25:54][,apply(data_hi3[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo3 <- paste("ISO_dunno0~", "-1 + as.numeric(own_foreign>0)+as.numeric(sales_export>0) +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo3[,25:54][,apply(data_lo3[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi3<-glm(formula_hi3, data=data_hi3, family=binomial(link=logit))
model_lo3<-glm(formula_lo3, data=data_lo3, family=binomial(link=logit))

#4. NIT split, country fe
data_hi4<-data[data$NITCORR>=median(tapply(data$NITCORR, data$COUNTRY, mean), na.rm=T),]
data_lo4<-data[data$NITCORR<median(tapply(data$NITCORR, data$COUNTRY, mean), na.rm=T),]

formula_hi4 <- paste("ISO_dunno0~", "-1 + as.numeric(own_foreign>0)+as.numeric(sales_export>0) +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi4[,25:54][,apply(data_hi4[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo4 <- paste("ISO_dunno0~", "-1 + as.numeric(own_foreign>0)+as.numeric(sales_export>0) +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo4[,25:54][,apply(data_lo4[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi4<-glm(formula_hi4, data=data_hi4, family=binomial(link=logit))
model_lo4<-glm(formula_lo4, data=data_lo4, family=binomial(link=logit))



#exports/foreign>0.1
#1. CORR08 split, country fe
data_hi1<-data[data$CORR08>=median(tapply(data$CORR08, data$COUNTRY, mean), na.rm=T),]
data_lo1<-data[data$CORR08<median(tapply(data$CORR08, data$COUNTRY, mean), na.rm=T),]

formula_hi1 <- paste("ISO_dunno0~", "-1 + as.numeric(own_foreign>0.1)+as.numeric(sales_export>0.1) +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi1[,25:54][,apply(data_hi1[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo1 <- paste("ISO_dunno0~", "-1 + as.numeric(own_foreign>0.1)+as.numeric(sales_export>0.1) +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo1[,25:54][,apply(data_lo1[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi1<-glm(formula_hi1, data=data_hi1, family=binomial(link=logit))
model_lo1<-glm(formula_lo1, data=data_lo1, family=binomial(link=logit))

#2. CPI split, country fe
data_hi2<-data[data$CPI08>=median(tapply(data$CPI08, data$COUNTRY, mean), na.rm=T),]
data_lo2<-data[data$CPI08<median(tapply(data$CPI08, data$COUNTRY, mean), na.rm=T),]

formula_hi2 <- paste("ISO_dunno0~", "-1 + as.numeric(own_foreign>0.1)+as.numeric(sales_export>0.1) +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi2[,25:54][,apply(data_hi2[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo2 <- paste("ISO_dunno0~", "-1 + as.numeric(own_foreign>0.1)+as.numeric(sales_export>0.1) +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo2[,25:54][,apply(data_lo2[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi2<-glm(formula_hi2, data=data_hi2, family=binomial(link=logit))
model_lo2<-glm(formula_lo2, data=data_lo2, family=binomial(link=logit))

#3. GOVEFF split, country fe
data_hi3<-data[data$GOVEFF08>=median(tapply(data$GOVEFF08, data$COUNTRY, mean), na.rm=T),]
data_lo3<-data[data$GOVEFF08<median(tapply(data$GOVEFF08, data$COUNTRY, mean), na.rm=T),]

formula_hi3 <- paste("ISO_dunno0~", "-1 + as.numeric(own_foreign>0.1)+as.numeric(sales_export>0.1) +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi3[,25:54][,apply(data_hi3[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo3 <- paste("ISO_dunno0~", "-1 + as.numeric(own_foreign>0.1)+as.numeric(sales_export>0.1) +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo3[,25:54][,apply(data_lo3[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi3<-glm(formula_hi3, data=data_hi3, family=binomial(link=logit))
model_lo3<-glm(formula_lo3, data=data_lo3, family=binomial(link=logit))

#4. NIT split, country fe
data_hi4<-data[data$NITCORR>=median(tapply(data$NITCORR, data$COUNTRY, mean), na.rm=T),]
data_lo4<-data[data$NITCORR<median(tapply(data$NITCORR, data$COUNTRY, mean), na.rm=T),]

formula_hi4 <- paste("ISO_dunno0~", "-1 + as.numeric(own_foreign>0.1)+as.numeric(sales_export>0.1) +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi4[,25:54][,apply(data_hi4[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo4 <- paste("ISO_dunno0~", "-1 + as.numeric(own_foreign>0.1)+as.numeric(sales_export>0.1) +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo4[,25:54][,apply(data_lo4[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi4<-glm(formula_hi4, data=data_hi4, family=binomial(link=logit))
model_lo4<-glm(formula_lo4, data=data_lo4, family=binomial(link=logit))




#####USE FOREIGN PRESSURE VARS

#1. CORR08 split, country fe
data_hi1<-data[data$CORR08>=median(tapply(data$CORR08, data$COUNTRY, mean), na.rm=T),]
data_lo1<-data[data$CORR08<median(tapply(data$CORR08, data$COUNTRY, mean), na.rm=T),]

formula_hi1 <- paste("ISO_dunno0~", "-1 + foreignpressure  +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi1[,25:54][,apply(data_hi1[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo1 <- paste("ISO_dunno0~", "-1 + foreignpressure  +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo1[,25:54][,apply(data_lo1[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi1<-glm(formula_hi1, data=data_hi1, family=binomial(link=logit))
model_lo1<-glm(formula_lo1, data=data_lo1, family=binomial(link=logit))

#2. CPI split, country fe
data_hi2<-data[data$CPI08>=median(tapply(data$CPI08, data$COUNTRY, mean), na.rm=T),]
data_lo2<-data[data$CPI08<median(tapply(data$CPI08, data$COUNTRY, mean), na.rm=T),]

formula_hi2 <- paste("ISO_dunno0~", "-1 + foreignpressure   +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi2[,25:54][,apply(data_hi2[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo2 <- paste("ISO_dunno0~", "-1 + foreignpressure   +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo2[,25:54][,apply(data_lo2[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi2<-glm(formula_hi2, data=data_hi2, family=binomial(link=logit))
model_lo2<-glm(formula_lo2, data=data_lo2, family=binomial(link=logit))

#3. GOVEFF split, country fe
data_hi3<-data[data$GOVEFF08>=median(tapply(data$GOVEFF08, data$COUNTRY, mean), na.rm=T),]
data_lo3<-data[data$GOVEFF08<median(tapply(data$GOVEFF08, data$COUNTRY, mean), na.rm=T),]

formula_hi3 <- paste("ISO_dunno0~", "-1 + foreignpressure   +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi3[,25:54][,apply(data_hi3[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo3 <- paste("ISO_dunno0~", "-1 + foreignpressure   +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo3[,25:54][,apply(data_lo3[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi3<-glm(formula_hi3, data=data_hi3, family=binomial(link=logit))
model_lo3<-glm(formula_lo3, data=data_lo3, family=binomial(link=logit))

#4. NIT split, country fe
data_hi4<-data[data$NITCORR>=median(tapply(data$NITCORR, data$COUNTRY, mean), na.rm=T),]
data_lo4<-data[data$NITCORR<median(tapply(data$NITCORR, data$COUNTRY, mean), na.rm=T),]

formula_hi4 <- paste("ISO_dunno0~", "-1 + foreignpressure   +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi4[,25:54][,apply(data_hi4[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo4 <- paste("ISO_dunno0~", "-1 + foreignpressure   +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo4[,25:54][,apply(data_lo4[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi4<-glm(formula_hi4, data=data_hi4, family=binomial(link=logit))
model_lo4<-glm(formula_lo4, data=data_lo4, family=binomial(link=logit))






#SIMULATION
#back to Model 1

data_hi1<-data[data$CORR08>=median(tapply(data$CORR08, data$COUNTRY, mean), na.rm=T),]
data_lo1<-data[data$CORR08<median(tapply(data$CORR08, data$COUNTRY, mean), na.rm=T),]

formula_hi1 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_hi1[,25:54][,apply(data_hi1[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")
formula_lo1 <- paste("ISO_dunno0~", "-1 + own_foreign+sales_export +logTIMEREG+how_old+firmsize+subsidiary+stateowned+SERVICE+",
	paste(names(data_lo1[,25:54][,apply(data_lo1[,25:54], 2,max,na.rm=T)==1]),collapse="+"),sep="")

model_hi1<-glm(formula_hi1, data=data_hi1, family=binomial(link=logit))
model_lo1<-glm(formula_lo1, data=data_lo1, family=binomial(link=logit))

sims <- 10000

pe_hi1<-coef(model_hi1)
vc_hi1<-vcov(model_hi1)
mdata_hi1 <- extractdata(formula_hi1, data_hi1, na.rm=TRUE)
simbetas_hi1 <- mvrnorm(sims, pe_hi1, vc_hi1)

pe_lo1<-coef(model_lo1)
vc_lo1<-vcov(model_lo1)
mdata_lo1 <- extractdata(formula_lo1, data_lo1, na.rm=TRUE)
simbetas_lo1 <- mvrnorm(sims, pe_lo1, vc_lo1)




#hi
xhyp1_hi1 <- cfMake(as.formula(formula_hi1), mdata_hi1, nscen=1, f=median)
xhyp1_hi1 <- cfChange(xhyp1_hi1, "own_foreign", x=1, xpre=0, scen=1)
yfd1_hi1 <- logitsimfd(xhyp1_hi1, simbetas_hi1)

xhyp2_hi1 <- cfMake(as.formula(formula_hi1), mdata_hi1, nscen=1, f=median)
xhyp2_hi1 <- cfChange(xhyp2_hi1, "sales_export", x=1, xpre=0, scen=1)
yfd2_hi1 <- logitsimfd(xhyp2_hi1, simbetas_hi1)


xhyp3_hi1 <- cfMake(as.formula(formula_hi1), mdata_hi1, nscen=1, f=median)
xhyp3_hi1 <- cfChange(xhyp3_hi1, "logTIMEREG",
x=mean(mdata_hi1$logTIMEREG)+sd(mdata_hi1$logTIMEREG), 
xpre=mean(mdata_hi1$logTIMEREG)-sd(mdata_hi1$logTIMEREG),
 scen=1)
yfd3_hi1 <- logitsimfd(xhyp3_hi1, simbetas_hi1)

results_hi1<-cbind(yfd1_hi1,yfd2_hi1,yfd3_hi1)
colnames(results_hi1)<-c("own_foreign", "sales_export",  "logTIMEREG")



#lo
xhyp1_lo1 <- cfMake(as.formula(formula_lo1), mdata_lo1, nscen=1, f=median)
xhyp1_lo1 <- cfChange(xhyp1_lo1, "own_foreign", x=1, xpre=0, scen=1)
yfd1_lo1 <- logitsimfd(xhyp1_lo1, simbetas_lo1)

xhyp2_lo1 <- cfMake(as.formula(formula_lo1), mdata_lo1, nscen=1, f=median)
xhyp2_lo1 <- cfChange(xhyp2_lo1, "sales_export", x=1, xpre=0, scen=1)
yfd2_lo1 <- logitsimfd(xhyp2_lo1, simbetas_lo1)


xhyp3_lo1 <- cfMake(as.formula(formula_lo1), mdata_lo1, nscen=1, f=median)
xhyp3_lo1 <- cfChange(xhyp3_lo1, "logTIMEREG",
x=mean(mdata_lo1$logTIMEREG)+sd(mdata_lo1$logTIMEREG), 
xpre=mean(mdata_lo1$logTIMEREG)-sd(mdata_lo1$logTIMEREG),
 scen=1)
yfd3_lo1 <- logitsimfd(xhyp3_lo1, simbetas_lo1)

results_lo1<-cbind(yfd1_lo1,yfd2_lo1,yfd3_lo1)
colnames(results_lo1)<-c("own_foreign", "sales_export",  "logTIMEREG")


varnames<-c("Foreign Ownership", "Sales for Export",  "Regulatory Burden")

trace_hi1 <- ropeladder(x=results_hi1[1,],
lower=results_hi1[2,],
upper=results_hi1[3,],
labels=varnames,
plot=1, lty=1, pch=19, 
sublabels=c("Strong Insts.","",""),
sublabelsyoffset=-0.05)

trace_lo1 <- ropeladder(x=results_lo1[1,],
lower=results_lo1[2,],
upper=results_lo1[3,],
labels=varnames,
plot=1, lty=3, pch=19, 
sublabels=c("Weak Insts.","",""),
sublabelsyoffset=-0.05)


vertmark <- linesTile(x = c(0,0),
y = c(0,1),
lty = "solid",
plot = 1
)

#create plot and save to PDF
tc <- tile(
	trace_hi1, trace_lo1, 
 vertmark,
output = list(width=6, file="isq_plot1", type="pdf"),   
xaxis=list(at=c(0,0.1,0.2,0.3), cex=0.9),
xaxistitle = list(labels="Change in Pr(Certification)"),
#maintitle = list(labels = "Substantive Effects in Different Institutional Contexts", cex=0.9),
gridlines=list(type="x", col="gray30")
) 



